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ABSTRACT 



in 



Context. Shortened version: 

Aims. The fate of interstellar clouds embedded in a hot tenuous medium depends on whether the clouds suffer from evaporation or 
C^h' whether material condensates onto them. Analytical solutions for the rate of evaporative mass loss from an isolated spherical cloud 

embedded in a hot tenuous gas are deduced by Cowie & McKee (1977). In order to test the validity of the analytical results for more 
^0 | realistic interstellar conditions the full hydrodynamical equations must be treated. 

Methods. Therefore, 2D numerical simulations of the evolution of interstellar clouds are performed with different internal density 
structures and surrounded by a hot plasma reservoir. Self-gravity, interstellar heating and cooling effects and heat conduction by 
electrons are added. Classical thermal conductivity of a fully ionized hydrogen plasma and saturated heat flux are considered. 
Results. Using pure hydrodynamics and classical heat flux we can reproduce the analytical results. Heat flux saturation reduces the 
evaporation rate by one order of magnitude below the analytical value, because the density distribution changes drastically during 
the simulation. This main result still holds with self gravity or different cloud density structures while keeping the cloud radius 
and surface temperature constant. The evolution changes, however, totally for more realistic conditions when interstellar heating 
and cooling effects stabilize the self-gravity. The evaporation then turns into condensation, because the additional energy by heat 
conduction can be transported away from the interface and radiated off efficiently from the cloud's inner parts. 
Conclusions. The consideration of a saturated heat flux is inevitable for interstellar clouds embedded in hot tenuous gas. Under more 
realistical conditions with radiative cooling heat conduction leads to condensation in contrast to analytical predictions which require 
evaporation. This provides an efficient way to accrete and mix intercloud material into clouds. 

^ ' Key words. ISM: clouds - ISM: structure - Conduction - Hydrodynamics - Method: numerics 
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• 1 . Introduction transfer. At this interfaces heat conduction plays an important 

• ' ,,»»,. /to«,n , , , ■ , role and is of prime importance whether the clouds evaporate or 

ON ■ The Interstellar Medium (ISM) can be described as an inh^mo- condense and therefore in mass Analytical approximations 

O , geneous ensemble of three phases (McKee & Ostnker \M&. ^ ^ ^ douds fa a hot ]m guffer from ration 

r- , The cold neutral phase with temperature T ~ 80 K and den- and ^ therefore fee Qn] shon M The ^ ^ interfjtel _ 

■ sity n ^ 40 cm is represented by the cores of molecular . , , pmhpA A p A in „ hot mpHi,™ arP „w™H ™ th» nth.r 



clouds which are confined by a warm neutral to slightly ionized hand ^ ^ m stab0ized inst evaporation so that 

medium (T ~ 8000 K, n ~ 0.3 cm 3 ). These two components ^ destraction is delayed 



are in pressure equilibrium if the gas is externally heated and „, , , , , , ,, ^ 

H ! can cool radiatively (Field et al. 1969). According to McKee ™ese observations include High-Velocity Clouds (HVCs) 

■ & Ostriker they are embedded in a third phase of the ISM: Hi entities consisting of a multi-phase structure (Wakker & 

, . j., . • . , , /TjTi/r\ -*u t A Schwarz 1991) and characterized by radial velocities that are 

the hot dilute intercloud medium (HIM) with T ~ 10° K and . '— ; — '. . . , , / , .... 

in-3 m-4 -3 c- tu- a ai incompatible with simple models of the differential rotation of 

n ~ 10 ° — 10 cm . Since this component is produced lo- , f . , , „, , , „ , . 

,, , , • , . . j the galactic disk (see Wakker & van Woerden ( 1997) tor a re- 

cally by supernova explosions at even higher temperatures and fe . . ^. . , * — r J 

, ... •. • • 11 t u • tu • vl.u cent review). Distance measurements of the cloud complexes 

densities, it can originally not be in pressure equilibrium with the . ..' , „ , „, ,. . _ F , . 

, ■ i -i jt c f. j i ., t-. ■ remain difficult. For at least two of them upper limits for their 

cooler phases and has therefore to expand vehemently. During ,. , , , , , , fT. t „ 

... • ■ , • , . , TTf A/ r t , . u u.u distances are deduced by Danly et al. (1993), Keenan et al. 

this expansion shocks arise and the HIM penetrates through the , , , rrrwr=n , . , , 

, • . , Tcnv/r x-, , , . , . , . (1995) and van Woerden et al. (1997) which locate them in 

ambient clumpy ISM. Denser clouds cannot be swept-up but , — , , . , , _ . . ' . • , , , 

a u Ji utiv/t ,i, t t u u u a a a t u ■ the hot galactic halo. Therefore, interactions with the hot rar- 

are passed by the HIM so that they become embedded therein. _ , , , 6 

n , . ■• c ., u ■ i . . . . tU ened halo gas are expected. Indications for such interactions are 

Due to the strong discrepancy of the physical states between the , . r , , • , • , 

, • . ? , . i , . , , . , the detection of morphological percuhanties like head-tail struc- 

phases an interface has to form where the hot phase and the i-WC™ F T . . , , , 

, , i j • t t t . a a tures (Bruns et al. 2000; 2001 ). Interactions with gas belonging 

molecular cloud are in contact. Temperatures and densities are , " — ; — : — — rf 1 — : — ' ... , , , ? 

, . .i u t a- t ,u , i a t a to the galactic disk and having lower velocities are deduced from 

connected through steep gradients that lead to energy and mass , , & . . „,,.,., jtt , 

the observation of so-called velocity bridges dPietz et al. 1996). 
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Further examples for interstellar clouds in a stream of hot gas 
are cometary globules, like G134. 6+1.4, which is located in a 
galactic chimney (Normande au et al. 19 96). This chimney is as- 
sociated with the Hll region W4 and forms a cone in the Hi 
layer above W4 with a diameter of about 100 pc. At the bot- 
tom of the chimney the starcluster IC 1805 with 9 O stars can be 
found ( |Heyer et al. 1996| |Tayloretal. 1999| l. An intense inter- 
action between the shocked hot stellar winds and G 134. 6+ 1.4 is 
expected. Energy exchange between the cold molecular to neu- 
tral gas of the globule and the hot chimney gas leads either to 
an evaporation and therefore a destruction of the cloud what is 
not observed or to condensation of hot gas onto the globule and 
therefore to a stabilization. 

Previous papers dealing with the evaporation and condensa- 
tion of molecular clouds only solve the time-independent en- 
ergy equation without taking temporal effects into account like 
mass changes, heating, etc., and dynamical effects like mixing 
of warm cloud material with the hot plasma and surface insta- 
bilities. In this paper we examine the evolution of interstellar 
clouds that are embedded in a hot, dilute medium by solving the 
full hydrodynamical equations numerically and furthermore fol- 
lowing the heat conduction process in detail. In §2 we discuss 
the predicted evolution that is expected refering to the papers of 
Cowie & McKee ( fl977l hereafter: CM77) and Dalton & Balbus 
d 19931 hereafter: DB93). The treatment of heat conduction in the 
context of hydrodynamical simulations is described in §3. In §4 
we introduce the different models and show the results of the 
simulations. Conclusions are drawn in §5. 



2. Heat Conduction 

The energy excange between the phases as a consequence of heat 
conduction is described by a heat flux q. In the classical case this 
flux is calculated assuming a diffusion appoximation leading to 
the formula given by Spitzer ( 1962)): 
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with the heat conduction coefficient 

1.84 x 10~ 5 T 5 / 2 . 
k = erg s K cm 

where the Coulomb logarithm is 



(1) 



(2) 



equation (|Khan & Rog nlien 1981} IMatte & Virmont 1982; 
ICampbell 1984f ~ 

The heat conduction coefficient k is a function of the ionization 
state of the medium. At temperatures lower than 10 4 K the ion- 
ization fraction of the gas is negligible and with this the heat 
conduction by electrons. Inside the cloud only conduction by 
neutrals may play a role for the temperature distribution. On the 
other hand the temperature gradient in this part of the cloud is too 
smooth to create a significant heat flux so that heat conduction 
becomes negligible. The evaporation and condensation of inter- 
stellar clouds was first investigated by Zel'dovich & Pikel'ner 
( 119691) and later by Penston & Brown ( 19701 1 in a plane parallel 
approximation. The more realistic case of spherical clouds was 
examined by Graham & Langer ( 1973 ). They solved the time- 
independent energy equation including heat conduction. As a 
main result they concluded that clouds with radii below a critical 
radius suffer from evaporation while clouds with larger radii gain 
mass by condensation of material onto the cloud surface because 
radiative cooling exceeds the heat input due to heat conduction. 
Numerical simulations were performed by Chevalier (119751 1 in 
order to describe the evolution of interstellar clouds in young 
supernova remnants. Although the temperatures rose up to about 
10 7 K he took only the classical heat flux into account and, by 
this, overestimated the mass-loss rate by some orders of magni- 
tude. 

CM77 made a general study of the evaporation of spherical 
clouds including saturation effects. In the treatment of the sat- 
uration of nonmagnetic, nonradiative, spherical evaporation, the 
domain outside the cloud is broken up into three zones, each of 
which is treated seperately. These zones are joined by prescribed 
boundary conditions for the heat flux and the temperature. The 
heat conduction of the innermost zone which inner boundary 
is the cloud surface is described by the classical heat flux. For 
larger radii the temperature rises and the heat flux becomes sat- 
urated. CM77 treated this by abruptly changing the mathemati- 
cal form of the conductivity at a well-defined saturation radius. 
The temperature rises steeply through the saturation zone un- 
til it reaches the value of the ambient gas. At this radius, CM77 
switched back to the classical description of the heat flux. Beside 
a local saturation parameter a as the ratio of the classical to the 
saturated heat flux, 
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with the electron density n e and the electron temperature Tg 
units of 10 6 K. 

This description breaks down if the local temperature scale- 
height falls below the mean free path of the conducting elec- 
trons. In this case the heat flux is replaced by a flux-limited form 
the so-called saturated heat flux (CM77): 



(4) 



with the sound speed c and density p. Q s is an efficiency 
factor less than or of the order of unity, which embodies 
some uncertainties connected with the flux-limited treatment 
and flux suppression due to magnetic fields. This method 
yields results in good agreement with laser-fusion experi- 
ments dMorse & Nielsen 1973] IMannheimer & Klein 19751 1 
and numerical simulations solving the Fokker-Planck 



they introduced a global saturation parameter ctq 



CTq = 



25$ s pc 3 i? 



T, 



1.54 x 10 7 K J n f $>,R 



1 



(6) 



in order to ascertain the influence of saturation effects only by 
knowing the cloud radius R or R pc the same in parsec and the 
heat conduction coefficient, the temperature and particle density 
of the ambient medium kj, Tf and rij. For <jq < 1 the mass-loss 
is given by 
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where k is Boltzmann's constant and p is the mean molecular 
weight of the ambient plasma. For larger values of cto the mass 
loss is a decreasing function of <tq (see Fig . [TJ . 
In a following paper McKee & Cowie i 19771 1 analysed the in- 
fluence of radiative cooling on the evaporation rate. They found 
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Fig. 1. Mass-loss rate m normalized to the classical value m c i ass 
(eq. as a function of the global saturation parameter ctq. 
Comparison of the rates refering to CM77 (dotted) and DB93 
(solid curve). 



out that the additional heat input due to heat conduction is com- 
pensated by radiative cooling for (To < 0.027/<I> s . Also in this 
study the authors devided the cloud environment into three zones 
where the physics are described in different ways. In the inner 
and the outer zone radiative losses are larger than the spherical 
divergence of the heat flux, while in the zone between the oppo- 
site is true and radiative losses are neglegted. Because of the fact 
that only the classical heat flux is considered, this approach can 
be only considered as a very crude estimation for the mass loss 
rate. Saturation reduces the energy input into the cloud so that 
radiation cooling can exceed this additional heat input. 
Begelman & McKee ( 119901) investigated the evolution of a two- 
phase medium consisting of clouds embedded in a hot plasma. 
They showed that heating and cooling processes dominate the 
hot phase and heat conduction can be neglected if the length 
scales of the relevant structures are larger than a critical lenght, 
the so-called Field length: 
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where Cm defines the maximum heating or cooling rate. The 
Field length is the maximum length scale on which energy trans- 
port by heat conduction is effective. So the temperature distribu- 
tion of structures with r ^ Xp that are embedded in the hot 
plasma is dominated by heat conduction while this process can 
be neglected for r > Af- For an isolated cloud they concluded 
that the classical mass-loss rate is valid for clouds with R <C Af 
while condensation occurs for R « 0.24 — 0.36Af. 
DB93 recovered the considerations of CM77 but this time used 
an effective heat flux that is the harmonic mean of the classical 
and the saturated flux in order to get a smooth transition between 
both regimes: 



dT 



<7eff = 



1 + er dr 



(9) 



Also in this description the time-independent energy equation 
can be solved analytically. This leads to the temperature distri- 
bution outside the cloud and also to the mass-loss rate as a func- 
tion of <jq. The mass loss rate is compared with the one derived 
by CM77 in Fig. [TJ The temperature distributions calculated by 




the formulae of DB93 for different values of Co (Fig. |2]i reflect 
the two different mathematical descriptions of the heat flux. The 
solutions for the classical flux (<7o < 1) show a steep tempera- 
ture increase at the edge of the cloud and a decreasing gradient 
for larger r. For large <7o the curvature is positive near the cloud 
surface so that Tj is reached at smaller r than in the classical 
case. 
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Fig. 2. Temperature profiles for different values of <tq derived 
fromDB93. 



Because only the energy equation is solved in CM77 as well 
as in DB93, dynamical effects like mixing of the cold cloud 
material with the hot plasma are not considered. On the other 
hand, it must be expected that the cloud material that is detached 
from the cloud surface is still cooler than the ambient plasma 
because only a slight rise in temperature of the cloud material 
leads to its overpressure that decouples the material from the 
cloud's surface. This will influence the temperature distribution 
near the surface and, by this, alter the mass-loss rate accordingly. 
Therefore, the hydrodynamical treatment is essential in order to 
explore to what extent the mass flow by evaporation and conden- 
sation and further gas-phase mixing changes the physical condi- 
tions around the interface. 



3. Hydrodynamical Treatment 

The spherical geometry of the cloud suggests to use a one di- 
mensional spherical symmetric hydro-code. With a limitation to 
only one spatial dimension it is impossible to follow the growth 
of instabilities that are expected when cold dense gas is accel- 
erated into warm rarefied gas. Therefore the evolution of clouds 
in a hot plasma is studied by two-dimensional hydrodynamical 
simulations. The Eulerian equations are solved on a rectangular 
cylindrically symmetric "staggered grid" which is of second or- 
der in space and based on the prescription of Rozyczka (1985). 
This code was extensively tested and used by different authors 
(e.g. Yorke & Welz 1996). The grid parameters, the resolution, 
the density /temperature profile of the initial cloud model and the 
considered physical processes are listed in Table[T]for the differ- 
ent models. 

The different setups for the homogeneous clouds are chosen in 
order to test the boundary conditions (model R2b) and the influ- 
ence of resolution on the mass-loss rate (model Rib, R2c). 
The boundary conditions at the upper, the left-hand and the right- 
hand sides are semi-permeable to allow for an outflow of gas 



Table 1. Model parameters and physical processes considered in the various simulations. For all models: i? c ioud = 41 pc, Tf = 
5.6 x 10 6 K, n f = 6.6 x 1CT 4 cm" 1 , a = 4.88. 
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class. 


class. & sat. 


gravity 


cooling 
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Rla 


800 x 400 
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homogen. 
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3.4 x 10 4 


Rib 


800 x 400 


0.5 


homogen. 
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3.4 x 10 4 


R2a 


200 x 100 
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homogen. 
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3.4 x 10 4 


R2b 


1000 x 500 
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homogen. 
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3.4 x 10 4 


R2c 


900 x 450 


0.33 


homogen. 




+ 






3.4 x 10 4 


R3 


200 x 100 
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core-halo 




+ 






6.4 x 10 4 


R4 


200 x 100 
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core-halo 




+ 


+ 




6.4 x 10 4 


R5 


200 x 100 


1 


core-halo 




+ 


+ 


+ 


6.4 x 10 4 



from the computational domain. The physical parameters at the 
lower boundary, the symmetry axis, are mirrored. In order to 
trace the condensation of the hot plasma onto the cloud a new 
quantity "colour" is introduced that is set in each cell to the den- 
sity fraction of hot ISM. At the beginning only the cells around 
the cloud possess a non-zero "colour" of value unity. During the 
calculation this quantity is advected like the others as e.g. mass 
or energy density. For the advection we use the monotonic trans- 
port (van Leer 1 1977b . The Poisson equation for self-gravity was 
solved for the corresponding model at each tenth timestep be- 
cause significant density changes of the cloud structure happen 
on a much larger timescale than the dynamical one. The energy 
equation includes heating, cooling and heat conduction, accord- 
ing to what physical process is included in the model: 



de 

dt 



V • (ev) = -PV v + T - A-V • q 



(10) 



Here e denotes the energy density, v the velocity, P the pressure, 
T the heating rate, A the cooling rate and q the heat flux. The 
equation of state for an ideal gas is assumed to be valid: 



P = (7 - l)e with 7 = 5/3 



(ID 



The used cooling function assumes collisional ionisation equi- 
librium and is a combination of the function introduced by 
Bohringer & Hensler (1989) for T > 1 4 K and solar metal- 
licity and by Dalgarno & McCray (119721 1 for the lower tempera- 
ture regime. The heating function considers cosmic rays (Black 
1987), X-rays and the photoelectric effect on dust grains (de 
Jong 119771 de Jong et al. [1980). This detailed description of 
heating and cooling processes cannot be incorporated in analyt- 
ical work like in CM77 that approximates the cooling function 
by four power-law segments. For temperatures above 10 4 K we 
used a more up-to-date description than CM77 who derived the 
cooling function in this temperature regime from Raymond, Cox 
& Smith (1976). Nevertheless, the cooling function used in our 
simulations is similar to the one used by CM77. 
The heat flux is calculated by taking both, the classical and the 
saturated flux into account. In order to apply a smooth transi- 
tion between classical and saturated regime we use the analytical 
form by Slavin & Cox ( fl992t 



Qsatl 1 ~ exp 
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This guarantees that the smaller flux is taken if both differ sig- 
nificantly. The heat flux due to electron diffusion is calculated 
separately using an implicit method which follows a scheme in- 
troduced by Crank & Nicolson (1947) and Juncosa & Young 



(I1971I) . To couple the two directions in space the method of frac- 
tional steps by Yanenko (119711 ) is used. A detailed description of 
the implementation as well as a numerical test of the heat con- 
duction code can be found in the online appendix. 
The simulated clouds are embedded in a hot plasma with a given 
temperature Tf = 5.6 x 10 6 K and particle density m = 6.6 x 
10~ 4 g cm~ 3 . The most realistic cloud model numbered as R5 
is generated for hydrostatic and thermal equilibrium under the 
constraint of spherical symmetry: 



pV ■ $ = -V • P 

r(r) = A(r) 



(13) 
(14) 



$(r) is the gravitational potential. The density and temperature 
profile of the cloud is then calculated by integrating equations 
(U~3b and (Tufl ) from inside-out using the core temperature of the 
cloud as a boundary condition and truncating the cloud's outer- 
most border where the energy density of the cloud reaches the 
value of the hot plasma ef. This condition is fulfilled for a ra- 
dius of i? c id = 41 pc. The density distribution we get in this 
way shows a clear core-halo structure with a tenuous rim and 
a pronounced central core (Figs. [3] and 2]) similar to observed 
interstellar clouds. 

This cloud radius together with the hot plasma parameters leads 
to a global saturation parameter of erg = 4.88 ($ s = 1) for 
which saturation of the heat flux is expected so that the classical 
mass-loss rate should be reduced for saturation effects according 
to CM77 by 51% to a value of 5.6 x 10 -4 M Q yr~ 4 or according 
to DB93 by 67% to a value of 3.7 x 1(T 4 M Q yr~\ respectively. 
In order to be close to the analytical considerations also a cloud 
with homogeneous density distribution is treated, denominated 
as Rl. With the same surface density and temperature as the 
stratified clouds it fulfills the condition of pressure equilibrium 
with the ambient medium. The heat conduction is here implied 
only classically. Self-gravity, heating or cooling are neglected. 
For models R2 the density and temperature structures are not al- 
tered but the physics are extended in a way, that now the effective 
heat flux is taken into account. These models can be used for a 
comparison with the analytical predictions of CM77 or DB93. 
For comparison reasons of all model clouds (also of the homo- 
geneous clouds) we fixed the radius to that of the R5 model, 
namely, 41 pc (Figs.[3]and[4]) and, by this, the value of co- 
in a more advanced model (R3) the density distribution is then 
modified according to a core-halo structure, while the temper- 
ature distribution has also to be changed in order to guarantee 
pressure equilibrium (Fig.|4](. The mathematical form of the heat 
conduction process remains the same. As the further improve- 
ment for model R4 we apply the same density structure as in 
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Fig. 3. Density profiles for the models with homogeneous den- 
sity distribution (dotted) and for those with pronounced central 
peak (solid line). 



model R3 but take now self-gravity into account. The tempera- 
ture profile has to change accordingly in order to serve for hydro- 
static equilibrium. And finally, as mentioned at the beginning of 
this part, R5 is developed as the most advanced model implying 
also thermal equilibrium. 



4.1. Model R1a-R1b 

Although the global saturation parameter cto suggests that satu- 
ration effects in the heat flux are to be expected, this first simu- 
lation is performed by only using the classical form of the heat 
flux. This is done in order to investigate the consequence of this 
simple approach on the mass-loss rate. The model starts with a 
homogeneous cloud of p = 8 x 10~ 24 g cm~ 3 and a spatially 
constant temperature of T = 1000 K. The temporal evolution 
of the density distribution is shown in Fig. [5] The velocity field 
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Fig. 4. Temperature profiles for the various models. Models 
Rl and R2 show a homogeneous distribution in order to fulfil 
pressure equilibrium with a homogeneous density distribution. 
Pressure equilibrium is necessary also for model R3 with a pro- 
nounced central density peak. The temperature profile for model 
R4 is calculated for hydrostatic equilibrium, model R5 has to 
fulfil thermal equilibrium as well. This does not alter the tem- 
perature distribution because of the chosen density profile. 



4. Models 

Here we present the evolution of the five generally different 
models and their variations with respect to the numerical res- 
olution. 



Fig. 5. Evolution of the density distribution of model Rla after 
Myr (a), 1 Myr (b), 2 Myr (c), 3 Myr (d), 4 Myr (e) and 5 Myr 
(f). Lines of equal density are drawn between p = 10~ 22 and 
10 -26 g cm -3 in steps of ldex. 

visualized by arrows in regions where p < 10 -25 g cm~ 3 is di- 
rected away from the cloud during the whole simulation. This 
is a clear indication for mass loss, while the inner parts of the 
cloud begin to collapse. The trigger for this process is the huge 
heat input at the cloud edge due to heat conduction, that rises the 
temperature there and leads to an overpressure with a two-fold 
consequence: First, it pushes material away from the surface into 
the ICM whereby the evaporated material is by far not as hot as 
the ICM when it decouples from the cloud. So the resulting tem- 
perature distribution outside the cloud is a rising function with 
the cloud distance (see Fig. |6j. Secondly, the evaporated mate- 
rial generates a repulsion that leads to a preassure wave running 
through the cloud and further to a pulsation of the cloud as a 
whole. The pulsation timescale corresponds very well to the dy- 
namical timescale of the cloud which is the sound crossing time. 
One has to keep in mind that the oscillations are excited artifi- 
cially by the switch-on of the numerical simulation. 
After some 10 5 years a quasistatic temperature distribution has 
formed (see Fig . |6j ■ In this figure also the pulsation of the cloud 
can be revealed. A comparison of T(r) at the end of the calcu- 
lation with an analytically derived distribution for ctq < 1 (see 
Fig. 13 shows a qualitatively good agreement. The steep rise in 
temperature at the cloud surface and a smooth asymptotic be- 
haviour for large r is typical for calculations with only the clas- 
sical heat flux. 

The mass-loss rate is calculated far away from the cloud at a 
distance of 95 pc using 

m = Airr 2 p(r)v{r) . (15) 




Fig. 6. Evolution of the temperature profile of model Rl. The 
course is in good agreement with an analytically derived distri- 
bution for <7o < 1 that takes only the classical heat flux into 
account. 



Positive values of rh stand for mass loss or evaporation of the 
cloud while negative values of rh indicate condensation of hot 
ISM onto the cloud. During the whole simulation evaporation 
occurs (see Fig. |7J and its strength is well correlated with the ra- 
dius of the cloud as expected from eq. [7] rh according to CM77 
or DB93 requires a value of about 1 x 10 -3 M yr _1 . This is in 
relatively good agreement with the value derived from the simu- 
lations (m w 0.8 x 10~ 3 M Q yr _1 ; see tabld2]). At the beginning 
rh represents the analytical value, and this more precisely for the 
spatially higher resolved simulation, what demonstrates the ac- 
curacy of our treatment. The model with higher resolution shows 




age [Myr] 

Fig. 7. Mass-loss rate for the models la (black solid curve) and 
lb (grey solid curve) compared with their cloud radius (lower 
black/grey dashed lines). 



no significant differences in its structure or evolution. Also the 
mass-loss rate is compareable. This is an indication that we can 
handle heat condution also in the low resolution case so that 



no time consuming highly resolved simulations are necessary 
as long as one is not interested in surface structural effects. 
The limitation of heat conduction to the classical case leads to 
a very large mass-loss rate in analytical studies as well as in 
numerical simulations. Due to the enormous energy input the 
evaporated material is accelerated to velocities near sonic speed. 

4.2. Model R2a-R2c 

These models differ from Rl in that saturation effects are now 
taken into account. The density and temperature structure of the 
initial model is the same as in Rl. The evolution in time of the 
density distribution is shown in Fig. [8] Because of the limitation 
of the energy transport due to the saturated heat flux, also the 
energy input and with this the rise in temperature at the cloud 
surface is reduced. The resulting overpressure at the cloud rim 
is lower than in model Rl. Nevertheless, a pressure wave that 
runs through the cloud and causes density fluctuations inside 
the cloud is triggered by the evaporated material (see innermost 
denser shell in Fig.[8j)). In contrast to model Rl, here no large- 
scale oscillations occur. 
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Fig. 8. Evolution of the density distribution of model R2a after 
Myr (a), 5 Myr (b), 10 Myr (c), 15 Myr (d), 20 Myr (e) and 25 
Myr (f). Lines of equal density are drawn between p = 10 
and 10 -26 g cm~ 3 in steps of 1 dex. 



22 



Also in this model cloud material is evaporated from the surface 
as a consequence of the additional heat input at the cloud surface 
due to heat conduction. Unlike model Rl the overpressure is not 
large enough to blow this material far and fast into space but it 
is mixed with the hot ICM and forms a transition layer around 
the cloud. The extension of this layer is visible in Fig. [H] by the 
growth of the distances between iso-density lines at the edge 
of the cloud. As time procedes Rayleigh-Taylor (RT) instability 
grows in the transition layer because warm and dense material is 
continuously accelerated into the hot rarefied plasma. The linear 

A where A is the 



growth timescale for RT instability is trt 
size of the initial perturbation which lies in the range of the cell 
size and a is the effective acceleration. Assuming that the evap- 
orated cloud material is accelerated to a value of 0.05 Mach over 
a distance of 3 pc in the transition layer, trt amounts to 0.8 
Myrs which is in good agreement with the value derived from 



the simulation. RT instability would be also expected to arise 
in models Rl, because the conditions for instability are analyti- 
cally fulfilled. However, the strong pulsations in this model are 
responsible to prevent its growth: phases in which RT instability 
is forced, i.e. when the cloud surface is accelerated outwards, are 
too short to allow its discernible growth. 

The extension of the transition layer is better visualized in Fig. [9] 
where the density profiles of the different models R2a-R2c are 
plotted as a function of the radius at time t = 30 Myr. The tran- 



CM77 cannot account for this physical process, the evaporation 
in their model remains at sound speed. The slight differences in 
the simulations with different spatial resolutions and grid exten- 
sion are not significant so that simulations with a resolution of 1 
pc cell - 1 are reliable. A comparison of the temperature profile of 
the simulation (see Fig. ITOb with an analytical one for Co = 4.8 
yields no good result whereas a theoretical profile for oo = 100 
fits the simulated model very well again. Conclusively, the tem- 
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Fig. 9. Density profiles of the models R2a-R2c at a time t = 
30 Myr. The transition region for 40 pc < r < 48 pc is clearly 
visible. 
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Fig. 10. Comparison of different temperature profiles between 
the model R2a after 62 Myr (solid curve) and theoretical profiles 
of constant 4.8 (dotted), = 100 (dashed). Additionally the 
local saturation parameter cr(r) deduced from the simulation is 
plotted (dashed-dotted). The temperature is normalized to Tf = 
5.5xl0 6 K. 



sition layer is clearly visible between 40 pc < r < 48 pc. The 
outer edge of the layer cannot be sharply defined. Also for larger 
radii up to 85 pc significant amounts of cloud material can be 
found. 

The results of the different R2 models show only weak depen- 
dencies on resolution or boundary conditions. In particular, the 
plateau between 40 pc and 45 pc in all R2 models in Fig. [^re- 
mains the same and reflects the regime where RT instability is 
excited. The fluctuations inside the cloud result from the pres- 
sure waves triggered by the repulsion of the evaporated mate- 
rial. The temperature inside the transition zone is rising with ra- 
dius and has values of T w 1.5 x 10 3 K at the inner boundary 
(r = 40 pc) and reaches T w 10 5 K at r =48 pc. Because 
the density structure of the layer runs oppositely to the temper- 
ature structure, the material in the transition zone can exist al- 
most in pressure equilibrium. Nevertheless, the remaining slight 
pressure gradient at the outer boundary leads to the observed 
evaporation. Like in the previous model the velocity field is di- 
rected outwards and indicates evaporation. In contrast to model 
Rl, the maximum speed of the lost material is subsonic with 
only 0.05 Mach (w 21 km s _1 ). This mass-loss rate amounts to 
rh ~ 1.4 x 10~ 5 M Q yr _1 (see Fig. [L3l and is therefore only 
1/40 of the analytically predicted one and agrees at best for a 
model with <tq w 100. We can show that the low Mach number 
is caused by the energy loss of the expanding (and evaporating) 
shell due to p ■ dV work against the external pressure. From our 
R2 model one can derive that the specific thermal energy is re- 
duced by a factor of almost 500. This value can be derived from 
the density and temperature values at the cloud edge (at 46 pc) 
and at r = 67 pc where Tt = 5.6 • 10 6 K is reached. The density 
contrast amounts to 3 • 10 3 (see Fig. [9]), the temperature con- 
trast to 0.16 (see Fig.fTOii. Since the analytical consideration by 



perature profile is the most important factor for the value of the 
evaporation rate and, on the other hand, results from the local 
heat flux that is a function of the density in the case of saturated 
heat conduction. With the formation of the transition zone at the 
edge of the cloud the density distribution has changed from the 
initial model and so does the evaporation rate. 

4.3. Model R3 

High-resolution observations of interstellar clouds, e.g. HVCs, 
show that most of them consist of an inhomogeneous gas and 
dust distribution with dense cores and rarefied envelopes. In or- 
der to analyze the influence of the heat conduction process on 
more realistic core-halo structures, the following cloud model 
consists of a pronounced core without self-gravity or heating 
and cooling. The temperature distribution of the initial model 
was chosen to account for pressure equilibrium. In the core re- 
gion with some 10 particles cm~ 3 the temperature drops to about 
10 K. Towards the edge it rises to T — 300 K at r = 10 pc and 
T = 1000 K at r — 30 pc. Because the heat flux is very small 
in the inner regions as a consequence of the low temperature, its 
influence will be limited to the cloud edge. The evolution of the 
density distribution is shown in Fig. [TT] Similar to model R2, a 
transtion region forms that is visible by the growth of the separa- 
tion of density isolines at the edge of the cloud. Also the growth 
of RT instabilies is discernible. 

In order to analyze the evolution at the cloud rim, a density pro- 
file at t — 37 Myr is plotted in Fig. [12] and compared with the 
initial model. Clearly visible is the transition region with similar 
expansion like model R2. This expansion excites sound waves 
travelling inwards, by this, dissipating energy, and rising the 
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Fig. 11. Evolution of the density distribution of model R3 after 
Myr (a), 5 Myr (b), 10 Myr (c), 15 Myr (d), 20 Myr (e) and 25 
Myr (f). Lines of equal density are drawn between p — 10 
and 10~ 26 g cm -3 in steps of ldex. 
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Fig. 12. Density profile of the model R3 at a time t = 37 Myr 
(solid curve) in comparison with the initial model (dashed). 



temperature by about 150 K in the region 30 pc < r < 40 pc 
whereas the core region remains nearly unaffected. The cloud's 
density distribution becomes quasistatic after the redistribution 
of the density so that a new pressure equilibrium is reached. The 
cloud structure is now comparable to model R2 except the fact, 
that the innermost part with r < 30 pc is almost inert and charac- 
terized by its huge density and low temperature. Also the density 
distribution inside the transition layer is comparable to the one 
of model R2 (see Fig. [9] and [12] for comparison). Because the 
evaporation rate is strongly determined by the processes in this 
transition layer, it is expected, that the mass-loss rate of both 
models should be the same (see Fig. [T3l . The reason for the 
slightly lower evaporation rate for model R3 can be found by 
the fact, that evaporated material can be provided only from the 
homogeneous part of the cloud within 30 pc < r < 40 pc while 
the innermost part of the cloud is inert against heat conduction. 
Conclusively no significant differences can be found in the evo- 
lution between a multi-phase cloud and a homogeneous cloud. 
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Fig. 13. Evaportion rates for the models R2a-R2c and R3. 



4.4. Model R4 

Clouds like those simulated in model R3 are massive enough so 
that self-gravity is no longer negligible. Model R4 was therefore 
calculated with self-gravity by solving the Poisson equation for 
the corresponding density distribution. The temperature of the 
initial model was set for hydrostatic equilibrium. The dense core 
has therefore to achieve a higher temperature than model R3 in 
order to prevent the cloud from gravitational collapse. The core 
temperature is about T = 230 K, rises up to T = 920 K at r = 
10 pc and reaches T = 1000 K at r = 30 pc. Heat conduction is 
very sensitive to a rise in temperature (k oc T 5 / 2 , see eq. [2J so 
that energy transport into the core region is expected. 
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The evolution of the density distribution is shown in Fig. [14] A 
transition region at the cloud edge has formed with an exten- 
sion almost constant after a formation time of about 15 Myr. 
The structure of this layer is much more regular than in the pre- 
vious models without self-gravity where RT instabilities domi- 



nate the region. Self-gravity is a process that stabilizes the cloud 
against RT instability (Murray et al. [1993): the effective accel- 
eration which the surface material experiences as the net effect 
of evaporation (like in models R2 and R3) vs. self-gravity points 
radially inwards as the density gradient does. 
Similar to model R3 the cloud can be devided into an inert core, a 
heat-conduction dominated zone with nearly homogeneous den- 
sity distribution and the transition zone. They are visible by the 
density profile after 53 Myr in comparison to the initial model 
(Fig. n~5l > and the corresponding temperature profiles (Fig. [T6l >. 
Because of gravitation the density profile of the intermediate 
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Fig. 15. Density profile of the model R4 at a time t — 53 Myr 
(solid curve) in comparison with the initial model (dashed). 
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Fig. 16. Temperature profile of the model R4 at a time t = 53 
Myr (solid curve) in comparison with the initial model (dashed). 

zone is not as fiat as the temperature profile but approaches hy- 
drostatic equilibrium. 

During the whole simulation energy is transported into the cloud 
at r > 18 pc due to heat conduction. This energy flow is con- 
nected with a mass transport because of the induced overpres- 
sure. This mass flow is neither locally nor temporally constant 
but corresponds to the temperature of the transition zone. To 
visualize this correlation, the mass flow and the temperature at 
r = 40 pc (inner boundary of the transition zone) are plotted in 



Fig- [El Similar to the previous models the evaporation of cloud 
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Fig. 17. Evolution of the mass-loss rate (solid curve) and the 
temperature (dotted) at the inner boundary of the transition zone 
(r = 40 pc) for model R4. 

material leads to a repulsion and therefore triggers a pressure 
wave inside the cloud. The consequence of the inward-travelling 
wave, namely, the negative m is visible in Fig. [17] (e.g.: t = 27 
Myr) and correlated with a rise in temperature of the transition 
zone. The wave interacts with the density discontinuity at r = 18 
pc (e.g.: t = 30 Myr), is reflected, travels outwards, reenters the 
transition zone and, with that, pushes cold cloud material into 
the transition zone (e.g.: t = 32.5 Myr). While material continu- 
ously flows into the transition layer the temperature there is fur- 
ther reduced (e.g.: 32.5 < t < 42 Myr). This transport is stopped 
as soon as a new hydrostatic equilibrium has been adjusted. The 
temperature fluctuations at the inner boundary of the transition 
zone also become visible at the temperature distribution outside 
the cloud by means of the evaporated mass flow. E.g. at r = 95 
pc the temperature varies by 2% correlated with the evaporation 
rate of about 1.5 x 10~ 4 M Q yr _1 on average (see Fig.[T8l. This 
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Fig. 18. Evaportion rate for the model R4. 

is slightly larger than in model R3 because the area feeding the 
evaporation is larger than in the previous model. 
This study shows that dynamical processes inside the cloud must 
not be neglegted but determine the transport inside the tran- 



sition zone and with this the temperature profile at the cloud 
edge. Even far away from the cloud this dynamics can be com- 
prehended by the mass-loss rate or temperature distribution. 
Because of the much larger sound velocity than that of the evap- 
oration flow the external temperature fluctuation affect the evap- 
oraton rate. 

4.5. Model R5 

Because already the studies by McKee & Cowie (1977) indicate 
that a reduction of the heat flux due to cooling processes is ex- 
pected but only for clouds with ctq < 0.027, model R5 takes 
heating and cooling processes into account. Opposite to CM77 
this model is calculated with self-gravity. From model R4 we 
conclude that self-gravity acts only as a smoothing agent for re- 
ducing instabilities that occur when heated material is evapo- 
rated from the cloud surface into the ICM. The mass-loss rate 
was not significantly altered by this process. So, this model 
should be comparable to the analytic model of Cowie & McKee. 
In models R2-R4 heat conduction leads to a temperature and 
therefore to a pressure increase at the edge of the cloud where the 
transition zone forms between the cold cloud and the hot ICM. 
In order to compensate the local energy input, two processes are 
physically plausible: The cloud distributes the energy through- 
out its volume by dynamical transport processes and/or it emits 
the energy by the process of radiative cooling. Cooling is very 
sensitive to the particle density and is most effective in regions 
with high density. The question is therefore whether the cooling 
efficiency is large enough to alter the formation or evolution of 
the transition zone and, with that, to change the evaporation rate. 
There are two main differences visible in the evolution of the 
density distribution of model R5 (Fig. fT9l> in comparison to the 
previous models. First of all, the velocity field outside the cloud 
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Fig. 19. Evolution of the density distribution of model R5 after 
Myr (a), 5 Myr (b), 10 Myr (c), 15 Myr (d), 20 Myr (e) and 25 
Myr (f). Lines of equal density are drawn between p = 10 -22 
and 10 -26 g cm -3 in steps on ldex. 

is directed towards the cloud that clearly indicates accretion of 
ICM onto the cloud's surface. Secondly, the density distribution 
inside the cloud remains nearly unaffected by heat conduction. 
Similar to the models R2-R4 a transition zone is formed but its 
extension is much narrower. This means that the influence of 



heat conduction is restricted to only a thin layer around the cloud 
with a thickness of only few parsecs. 
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Fig. 20. Comparison of temperature profiles of model R5 at dif- 
ferent ages: initial model (solid curve), after 30 Myr (dotted), 
after 40 Myr (dashed) and after 55 Myr (dashed-dotted). The 
radii (r*) at which the energy input due to heat conduction plus 
radiative heating is balanced by radiative cooling are marked by 
crosses. 



The used cooling function shows a steep increase for T around 
10 4 K. For regions with temperatures below 10 4 K this means 
that they can be heated up to about 10 4 K with high efficiency 
but then should somewhere easily find an energetic balance be- 
tween heating and cooling. This equilibrium situation is deter- 
mined by the Field length Ap (Eq. [8]). In the outermost shell 
of the cloud Ap amounts to 1.4 x 10 3 i? c id (Vieser & Hensler 
2005; hereafter: Paper II). For this reason the occurence of RT 
instability is prevented and no further refined grid resolution is 
neccesary. Due to the density gradient towards the cloud core 
Ap drops within an interface formed by heat conduction down to 
values much smaller than i? c id- Because of its high cooling effi- 
ciency by larger density a temperature of 8000- 10 4 Kis achieved 
at which its thermal energy density is not exceeding the gravi- 
tational one. As a consequence, the interface does not dissolve 
from the cloud but accumulates surrounding gas that enters the 
interface by means of heat conduction. 

From Fig. |20]one can recognize that between 30 Myrs and 40 
Myrs a quasi-stationary temperature distribution has been ad- 
justed. The balance at which the energy input due to heat con- 
duction plus heating is compensated by cooling is reached at 
temperatures of 10 4 K and radii r*. For the reached values at 
r+ the Field length amounts to not more than 10~ 2 i? c u- This 
effect has two consequences: At first, the physical one is that 
the parameter regime of condensation is reached; secondly, the 
cell size is larger than the typical thermodynamical length scale, 
i.e. here the mean free path of electrons and ions, respectively. 
Therefore, this means that the hydrodynamical treatment is ap- 
propriate (Koyama & Inutsuka 2004). 

Thermal energy is transported from r > inwards during the 
whole calculation and finds an equilibrium with the cooling rate 
at r* « 41 pc and 8000 K. Consequently, condensation occurs 
(Fig. [2Tb in this model with a value of about —3.5 x 10~ 6 M Q 
yr . Although the available conductive energy at r < r* is 
reduced the temperature rises moderately above 10 3 K for r > 
38 pc (Fig. [20b. by this generating a turbulent velocity field that 








"8 I i i i i I i i i i I i i i i I i i i i I i i i i L_ 

35 40 45 50 55 

age [Myr] 

Fig. 21. Evolution of the mass-loss rate of model R5. 
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Fig. 22. Evolution of the ICM mass fraction (/?ism/ Pc\d) of 
model R5 after Myr (a), 10 Myr (b), 20 Myr (c), 30 Myr (d), 
40 Myr (e) and 50 Myr (f). The solid line represents the cloud 
boundary. 



transports the accreted material from the cloud surface towards 
deeper layers. Due to the condensation of ICM onto the cloud 
and its transport inside the cloud, the mass fraction of ICM is 
rising continuously. The distribution of ICM in comparison to 
the whole cloud mass is visualized in Fig. [22] A growing shell 
becomes visible with an ISM mass fraction of at least 1%. After 
50 Myr this layer has reached a thickness of about 10 pc. 
To emphasise, heat conduction together with cooling processes 
does not only change the mass flow from evaporation to conden- 
sation but also provide a new possibility to enrich the cloud with 
ICM and to distribute it throughout large parts of the cloud. 

5. Discussion 

In this paper we investigate the evolution of interstellar clouds 
that are resting in a hot rarefied medium. Models with different 
density structures and various physical processes are compared 
with the aid of numerical simulations. The classical evaporation 
rate without taking saturation effects into account amounts to 



10 3 M Q yr 1 . Analytical consideration of saturation allows to 
lower it to about 5.6 x 10" 4 M yr" 1 (CM77) or 3.7 x 10" 4 
M Q yr -1 (DB93), respectively. 

Table 2. Mass-loss rates of the various simulations (negative val- 
ues mean mass accretion!) 
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The global saturation parameter of heat conduction in all calcu- 
lations here was <tq — 4.88. The mass-loss rates deduced from 
the simulations are summarized in Table [2] In order to com- 
pare with analytical studies (CM77, DB93), the initial model Rl 
shows a homogeneous density distribution and a temperature de- 
termined by pressure equilibrium. Influences by self-gravity or 
heating and cooling are neglected. As a first step, in model Rl 
only the classical description of heat conduction is taken into 
account. The mass-loss rate is in good agreement with the ana- 
lytical results, but we learn that it oscillates due to excited cloud 
pulsations. 

The evolution of a cloud taking saturation effects into account 
(model R2) proceeds not as violent as in Rl. The energy in- 
put into the cloud is reduced due to the saturated heat flux. 
Nevertheless the cloud's outermost layer is heated up and forms 
a transition layer between the cold cloud and the hot ICM in 
which the material can remain near pressure equilibrium with 
the ICM. From there cloud gas is dissolved and pushed into the 
hot ICM at subsonic speed. This, however, changes the temper- 
ature and density profiles in the cloud's environment drastically 
in comparison to the initial model, so that it is not surprising that 
the mass-loss rate is reduced to about one order of magnitude 
below the analytical value (CM77, DB93). On the other hand, 
this value is similar to the one derived analytically for do « 100 
and indeed the mean local saturation parameter deduced from 
the simulated temperature profile yields about 100 in contrast to 
4.88 set as the formal value. From that we learn, that the global 
saturation parameter is therefore not conclusive to derive the 
mass-loss rate. 

The influence of a multi-phase density structure of the cloud on 
the evaporation rate is not significant (model R3). Since density 
and temperature at the cloud surface are similar to the previous 
model which serve as the dominant factors for the heat input and 
by this for the mass loss rate, the same result is expected. 
That RT instability is the responsible mechanism that shapes the 
clouds' surfaces in models R2 and R3 becomes plausible, be- 
cause in model R4 the acceleration by self-gravity overcomes 
the outwards directed acceleration of the expanding shell. Thus, 
the shape of the cloud remains undisturbed but the mass-loss rate 
is slightly higher than in model R3 because the temperature near 
the cloud core of the initial model is higher for hydrostatic equi- 
librium. By this, heat conduction can affect deeper layers of the 
cloud. 

In model R5 heating and cooling processes are added. Since 
the used cooling function increases steeply at 10 4 K the heat- 
ing due to heat conduction can be balanced by radiative cooling 
at around 10 4 K for sufficiently high densities. In the presented 
cloud model saturated heat conduction is not able to transport 
enough thermal energy into the cloud in order to rise the temper- 



ture above 10 4 K and therefore to evaporate the cloud. The cloud 
rim remains at temperatures of about 10 4 K. Correlated with the 
energy transport onto the cloud's surface, mass is accumulated 
from the ISM and condenses onto the cloud. The cloud regions 
near the surface are slightly heated so that a turbulent velocity 
field emerges that mixes the accreted ICM into deeper layers. 
Three major conclusions can be drawn from the numerical mod- 
els presented in this paper: 1) Evaporation itself alters the envi- 
ronmental conditions in the sense of positive feedback. 2) The 
assumption of classical heat conduction is invalid under real 
conditions of interstellar clouds. In contradiction to those ana- 
lytical results, heat conduction is reduced to the saturated limit. 
3) In connection with radiative cooling this can lead to conden- 
sation of hot ICM onto clouds in a regime of cloud parameters 
where evaporation is requested from the analytical approach. By 
this it provides an alternative way to accrete and mix ICM with 
cloudy material. 
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Online Material 



Appendix A: Heat conduction - implementation & 
tests 

The implementation of heat conduction into the existing and 
tested hydro-code turned out to be more complicated than ex- 
pected in order not to limit the hydro-timestep Th y d to the 
conduction-timestep given by 



g 2 

Tcond = — [min(dr, dz)} 



Tk 



(A.l) 



where dr and dz are the distances between the centers of two 
neighbouring cells in r and z direction respectively. The problem 
dealing with these incompatible timescales is treated by using a 
semi-implicit scheme to solve the heat conduction equation. 
The equation describing the change of the internal energy den- 
sity due to heat conduction 



de 
dt 
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r dr 
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(A.2) 



is transformed via the equation of state for an ideal gas into an 
expression for the change of the temperature: 



dT T 1 d 
dt e r dr 
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dr 
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e dz 
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(A3) 



Both spatial directions are treated seperately using the same 
strategy. Therefore only the solution of the radial part is de- 
scribed in the following: The new temperature distribution is cal- 
culated for all cells in the r-direction with a fixed z-coordinate. 
This is done using an ADI-method. After the temperature in this 
row has converged, the index of the z-coordinate is increased by 
one and the equation is solved for the next row. 
If the discretisation of eq. IA.3I is done in a pure explicit 
way not only the timestep is limited to the conduction- 
timestep. This procedure by itself is numerically unstable 
(Crank & Nicolson I19471 I. Instead of this approach we used a 
semi-implicit procedure where the implicit part is weighted by a 
factor 0.5 < a < 1 and the explicit part by a factor 1 — a what 
results in: 



rpn 



A* 



(1 



2 TJj J_ 

ry+ry+i eA Ar 

Till rpn^ 



Ar 



rp II rp n 



A? 



/Tin+l 



•fij+i e^ +1 Ar 



(A.4) 



/ T-in+l rp 

K i,j-t-l 



y 




This equation can be rearranged and the different coefficients be 
combined to form: 



(A.5) 



For a row with the index i kept constant this describes a max- 
trix equation which can be solved with standard methods (e.g. 



Press et al. 119861 1. We used an ADI-method. The weightning 
parameter a is chosen in a manner so that the implicit part is 
pronounced if the conduction-timestep is much smaller than the 
hydro-timestep. This means that high temperatures or large con- 
ductivities dominate the integration domain. The ratio between 
the two timesteps is defined as (3 



r hyd 



that deals as a quantity to derive a: 

_ l-0-exp[-j3] 
a ~ /3(exp[-/3]-l) ' 

with the additional limitation 

a < 0.5 
0.5 < a < 1 
a > 1 



(A.6) 



(A.7) 



(A.8) 



so that the minimal implicit portion is 50%. The course of a as 
a function of (3 is plotted in Fig. IA.ll 



1.00 



0.90 r 



0.80 r 



0.70 r 



0.60 r 



0.50 r 




100.0 



Fig. A.l. Weightning factor a as a function of (3. 

The temperature distribution calculated is this way results from 
the assumption that the conductivities remains constant during 
the iteration of eq. IA.4I This is not the case as k is strongly tem- 
perature dependent. Therefore, after the new temperature distri- 
bution is calculated, the conduction coefficients have to be up- 
dated for this new temperatures and eq. IA.4I has to be solved 
again. This procedure is repeated until the error between the new 
temperatures and the old ones is less then 10~ 4 . 
We have tested the method described above through the compu- 
tation of problems with known solutions. Two test cases were 
run, one with a constant k and one with n cx T m with m = 5/2 
for which a self similar solution exists (Zel'dovich et al. 1969 ). 
The quality of the results for both test cases are similar so only 
the results of the k =const.-case are presented here. 
In this test the plasma is assumed permanently static and with a 
uniform and "frozen" density distribution with a particle density 
n = 6.6 x 10~ 4 cm~ 3 . We have considered the propagation of 
a spherical pure conduction front in a uniform medium. For this 
problem an analytic solution for the temperature distribution as 
function of time is available: 



T(r,t) = 



Q 



Airnt 



1.5 



exp 



4Qi 



(A.9) 



We have calculated the numerical solution in a cylindrical coor- 
dinate system, taking as initial profile the analytical solution at 
t = 0.5 s, for Q = 10 54 Kcm~ 3 andK = 1.36x 10 12 whatcorre- 
sponds to a conductivity of a plasma at a temperature of 5.6 x 10 6 
K. The grid is 200 x 100 cells with -0.032 < z < 0.032 pc and 
< r < 0.032 pc. The resulting propagation, compared with the 
analytical solution, and the relative error is shown in Fig. lA.2l for 
different times. The relative error within the hottest region of the 
plasma amounts to about 2-3%. Only at the edge of the conduc- 
tion front where the temperature is some orders of magnitudes 
lower than in the inner regions the error is slightly increased. In 
Fig. IA.2ti -f the influence of the grid boundary is visible. Here 
the temperature gradient is forced to be zero. With these studies 
we could verify the validity and accuracy of the code. 




Fig. A.2. Propagation of the conduction front after 1.0s (a), 6.0s (b), lis (c), 16s (d), 21s (e) and 26s (f). Comparison of the numerical 
results (diamonds) with the analytical solution (solid line). Additionally, the relative error between both results is plotted (convex 
line). 
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